\chapter{方程求根}
方程求根，即求解 
\begin{equation}
f(x)=0,
\end{equation}
实际计算中经常遇到没有解析解的方程，只能求数值解，比如
\begin{eqnarray}
\cos x -x =0,
\end{eqnarray}
只能求数值解。
有时$f(x)$甚至没有解析形式，可能是一系列计算的结果。

\section{二分法}	
已知区间 $(a,b)$ 内有 $f(x)=0$ 的根，且 $f(a)*f(b)<0$，可以如下逼近方程的根。
\begin{itemize}
\item [1] 取$c=(a+b)/2$，计算$f(c)$。
\item [2] 若$f(c)f(a) \leq 0$，说明$[a,c]$ 区间内存在根。将$c$的值赋给$b$，回到步骤１。
\item [3] 否则，即$f(c)f(a) > 0$，说明$f(c)*f(b)<0$，区间$[c,b]$内存在根，将$c$的值赋给$a$，回到步骤1。
\end{itemize}
每次迭代，区间都缩小一半。如此循环，不断缩小“嫌疑区间”，就可以逼近根的真实值。
当区间小于目标精度时，即可取区间的中点，作为方程的根。
%由于每次循环中区间的长度是上一次循环的一半，所以收敛速度是线性的。
我们可以写一个代码，实现这个算法
\begin{lstlisting}
#include<iostream>
#include<cmath>
#include<fstream>
#include<string>
using namespace std;

#define max_iteration 1E6

/*
 * findroot: 在[a,b]区间内寻找 func(x)=0　的根，并返回根的值
 * 算法：二分法
 */
double bisection(double (*func)(double), double a, double b, double precision){

        double ya=func(a), yb=func(b);
        double c=(a+b)/2, yc;
        int i;

        if( ya* yb > 0){//如果 func(a) * func(b) >0, 则[a,b]区间内不一定有解，拒绝任务
                cout<<"error: invalid input [a,b] for findroot: f(a)*f(b)>0 "<<endl;
        }

        for(i=0;i<max_iteration;i++){

                c = (a+b)/2;//区间中点
                yc = func(c);

                if( yc * ya <= 0 ){//如果 yc*ya<=0，说明[a,c]中有一个根
                        b=c;
                }
                else{//如果 yc*ya>0，说明[c,b]中有一个根
                        a=c;
                }
                if( fabs(a-b) < precision ){
                        return (a+b)/2;
                }
        }
        cout<<"bisection method: failed to find a root, after "<<max_iteration<<" iterations.\n";
        exit(1);
}

double g(double x){

        return x*x - 4*x -5;
}

int main(){

        double a=0, b=1, c;

        while( g(a) * g(b) >0 ){
                cout<<"Input valid boundaries inside which we get a root."<<endl;
                cout<<"a=";
                cin>>a;
                cout<<"b=";
                cin>>b;
        }

        cout<<" The root is x="<<(c=bisection(g,a,b,1E-9))<<" f(x)="<<g(c)<<endl;
        return 0;
}
\end{lstlisting}
程序中$g(x)$的形式可以任意改。
如果初始区间无效，程序会追问，直到得到有效的初始区间$[a,b]$。

\section{牛顿法}
如果方程的根是$s$，$f(s)$可以对任意$x$处做泰勒展开
\begin{equation}
0 = f(s) = f(x) + f'(x)(s-x) + \frac{ f''(x) }{ 2 } (s-x)^2 + \cdots
\end{equation}
当$s$很接近$x$时，忽略高阶项，有
\begin{equation}
s \approx x - f(x)/f'(x),
\end{equation}
受这个启发，可以构造一个迭代式
\begin{eqnarray}
x_{n+1} = x_n - \frac{ f(x_n) }{ f' (x_n)}
\end{eqnarray}

这个迭代具有几何的意义。
在$y=f(x)$曲线上$(x_n, f(x_n))$处作切线，与$x$轴的交点正是上面定义的$x_{n+1}$。
所以牛顿饭又叫做切线法。
迭代下去，当$f(x_{n+1})$小于某个很小的数时，即可把$x_{n+1}$当做方程根的近似值。
\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{NewtonMethodDemon.eps}
\caption{牛顿法逼近方程的根：切线与x轴交点为下一个迭代值。}
\label{figure:NewtonMethodDemon}
\end{figure}

\paragraph{收敛速度}
若真实值为 $\alpha$，$x_n$ 的误差为 $\epsilon_n = x_n - \alpha$，则
\begin{eqnarray}
x_{n+1} &=& x_n - \frac{ f(x_n) }{ f' (x_n)} \nonumber\\
				  &\approx& \alpha + \epsilon_n - \frac{ \epsilon_n f'(\alpha) + \frac{1}{2} \epsilon^2_n f''(\alpha)}{f'(\alpha) + \epsilon_n f''(\alpha)} \nonumber\\
				  	&\approx& \alpha + \frac{1}{2} \frac{ f''(\alpha) }{ f'(\alpha)} \epsilon^2_n,
\end{eqnarray}
即
\begin{eqnarray}
\epsilon_{n+1} \approx \frac{1}{2} \frac{ f''(\alpha) }{ f'(\alpha)} \epsilon^2_n,
\end{eqnarray}
所以，如果$\epsilon_n$是$10^{-3}$，$\frac{1}{2} \frac{ f''(\alpha) }{ f'(\alpha)} \approx 1$的话，$\epsilon_{n+1}$约为$10^{-6}$量级。
一般来说，牛顿法收敛得比二分法快。

但是牛顿法不一定收敛，一个典型的例子是，$f'(x_n) = 0$，迭代公式就没法用了。
可以先用二分法找到一个比较小的区间，再用牛顿法求方程的根。

另外，牛顿法需要使用一阶导数值，如果$f'(x)$没有解析值，而需要通过差分，可能引入新的误差，这种情况下，最好不要用牛顿法，可以使用后面介绍的割线法。

\begin{lstlisting}
/*
 * Newton's method: given a function (*func), which gives value and also the 1st derivative,
 * a suspicious zone [a,b], it finds a root in the zone within precision, and returns the root.
 */
double Newton_findroot(void (*func)(double, double &, double &), double a, double b, double precision){

        double x, dx, y, dy;
        int i;

        x = (a+b)/2;
        
        for(i=0;i<max_iteration;i++){
                func(x, y, dy);
                if( fabs(dy)<1E-9 ){
                        cout<<" error: f'=0 in Newton's root-finding method. \n";
                        exit(1);
                }
                dx = y/dy;
                x -= dx;
                if( (x-a)*(x-b) >0 ){
                        cout<<" error: x gets out of the given zone [a,b] in Newton's method.\n";
                        exit(1);
                }
                if( fabs(dx) < precision ){
                        return x;
                }
        }
        cout<<" Newton's method: after "<<max_iteration<<" steps, failed to find a root.\n";
        exit(1);
}
\end{lstlisting}

\section{割线法}
如果$f'(x)$很难得到，则可以用割线代替牛顿法中的切线。
\begin{itemize}
\item [1] 取两个初始点$x_0, x_1$，通过$(x_0, f(x_0)), (x_1, f(x_1))$两点作直线，交于 $x$ 轴，记该点为 $x_2$。
\item [2] 取点$x_1, x_2$，通过$(x_1, f(x_1)), (x_2, f(x_2))$两点作直线，交于 $x$ 轴，记该点为 $x_3$。
\item [3] $\cdots$
\end{itemize}
迭代公式为
\begin{eqnarray}
x_{n+1} = x_n - \frac{x_n - x_{n-1}}{f(x_n) - f(x_{n-1})} f(x_n).
\end{eqnarray}
\begin{figure}[h]
\centering
\includegraphics[width=0.7\textwidth]{SecantMethodDemon.eps}
\caption{割线法逼近方程的根：割线与x轴交点为下一个迭代值。}
\label{figure:SecantMethodDemon}
\end{figure}
可以用如下代码实现
\begin{lstlisting}
#include<iostream>
using namespace std;
#include<cmath>

#define max_iteration 1E6

double g(double x){
        return x*x-1;
}

/*
 * SecantMethod: given function func, suspicious zone [a,b], returns a root of func(x)==0
 */
double SecantMethod(double (*func)(double), double a, double b, double precision ){

        double c, fc, d, fd, e, fe;
        int i;

        c = (a+b)/2;//set starting points
        fc = func(c);
        d = (c+b)/2;
        fd = func(d);

        for(i=0;i<max_iteration;i++){
                e = d - (d-c)/(fd-fc)*fd;
                if( (e-a)*(e-b) >0 ){
                        cout<<"SecantMethod: error! iteration gets out of the given zone [a,b].\n";
                        exit(1);
                }
                fe = func(e);
                if( fabs(e-d)<precision && fabs(fe)<precision ){
                        return e;
                }
                c=d;
                fc=fd;
                d=e;
                fd=fe;
        }

        cout<<"After "<<max_iteration<<" iterations, failed to find one root"<<endl;
        exit(1);
}

int main(){

        double a=0, b=10;
        cout<<" root: x="<<SecantMethod(g, a, b, 1E-6)<<endl;
        return 0;
}
\end{lstlisting}

\paragraph{收敛速度}
记方程的根的真实值为 $\alpha$，$x_{n-1}, x_n$ 的误差为 $\epsilon_{n-1} = x_{n-1} - \alpha, \epsilon_n = x_n - \alpha$，则
\begin{eqnarray}
x_{n+1} &=& \alpha + \epsilon_n - \frac{ \epsilon_n - \epsilon_{n-1} } { f'(\alpha) (\epsilon_n - \epsilon_{n-1}) + \frac{ f''(\alpha) }{2} (\epsilon^2_n - \epsilon^2_{n-1})}
\nonumber\\
&\approx& \alpha + \frac{ f''(\alpha) }{2f'(\alpha)} \epsilon_n \epsilon_{n-1},
\end{eqnarray}
即
\begin{equation}
\epsilon_{n+1} \approx \frac{ f''(\alpha) }{ 2f'(\alpha) } \epsilon_n \epsilon_{n-1}.
\end{equation}
假设$ \epsilon_1 = \epsilon_0^r $，则
\begin{eqnarray}
\epsilon_0 &=& \epsilon_0^{1} \nonumber\\
\epsilon_1 &\propto& \epsilon_0^{r} \nonumber\\
\epsilon_2 &\propto& \epsilon_0^{r +1} \nonumber\\
\epsilon_3 &\propto& \epsilon_0^{2r +1} \nonumber\\
\epsilon_4 &\propto& \epsilon_0^{3r +2} \nonumber\\
\epsilon_5 &\propto& \epsilon_0^{5r +3} \nonumber\\
\epsilon_6 &\propto& \epsilon_0^{8r +5} \nonumber\\
\epsilon_7 &\propto& \epsilon_0^{13r +1} \nonumber\\
\cdots
\end{eqnarray}
其中出现了两列裴波那契数列 $1, 1, 2, 3, 5, 8, \cdots$，
显然有
\begin{eqnarray}
\epsilon_n \propto \epsilon_0^{a_{n-1} r + a_{n-2}}.
\end{eqnarray}
而裴波那契数列的通项公式为 
\begin{eqnarray}
a_n = \frac{1}{\sqrt{5}} \left[ (\frac{1+\sqrt{5}}{2})^n - (\frac{1-\sqrt{5}}{2})^n \right],
\end{eqnarray}
且
\begin{eqnarray}
\lim\limits_{n \rightarrow \infty} \frac{ a_n} {a_{n-1}} = \frac{ 1+ \sqrt{5}}{2} \approx 1.618033,
\end{eqnarray}
所以 $n \rightarrow \infty$ 时，
\begin{eqnarray}
\epsilon_{n+1} \propto \epsilon_{n}^{1.618033}.
\end{eqnarray}
显然，割线法没有牛顿法收敛得快。

\section{Stewenson方法}
显然，这个套路可以稍微再设计以下，比如Stewenson方法，
\begin{eqnarray}
x_{n+1} = x_n - \frac{ f(x_n)^2 }{f(x_n) - f(x_n - f(x_n))}.
\end{eqnarray}
即用 $(x_n, f(x_n)), (x_n - f(x_n), f(x_n - f(x_n)))$ 两点作割线，与割线法类似。
但是 Stewenson 方法收敛更快，下面我们给出证明。
假设方程的根的精确值为 $\alpha$，$x_n$ 的误差为 $ \epsilon_n = x_n - \alpha$，
经过推导，可以证明
%\begin{eqnarray}
%&& x_{n+1} \approx \alpha + \epsilon_n \nonumber\\
%&& - \frac{ (f'(\alpha) \epsilon_n + \frac{1}{2} f''(\alpha) \epsilon^2_n)^2 }
%{ f'(\alpha) \epsilon_n + \frac{1}{2}f''(\alpha)\epsilon^2_n - \left[ f'(\alpha)(\epsilon_n -f'(\alpha)\epsilon_n -\frac{1}{2}f''(\alpha)\epsilon^2_n) 																																	+ \frac{1}{2}f''(\alpha)(\epsilon_n - f'(\alpha)\epsilon_n -\frac{1}{2}f''(\alpha)\epsilon^2_n)^2 \right]}.																																	\nonumber\\
%\end{eqnarray}
%上式最后一项分子取到 $\epsilon^3_n$ 项，分母取到 $\epsilon^2_n$ 项，得到
%\begin{eqnarray}
%x_{n+1}  \approx \alpha + \frac{ f''(\alpha) (1- f'(\alpha)^2)}{f'(\alpha)^2} \epsilon^2_n,
%\end{eqnarray}
%即
%
\begin{eqnarray}
\epsilon_{n+1} \approx  \frac{ f''(\alpha) (1- f'(\alpha)^2)}{f'(\alpha)^2} \epsilon^2_n.
\end{eqnarray}
即Stewenson方法是二阶收敛的。
所以它不需要用到导数，但也可以和牛顿法收敛得差不多快。
缺点是公式稍微复杂一点点。

\section{勒让德多项式的根}
\subsection{勒让德多项式}
勒让德多项式可以通过递推式求得\cite{Arfken2005Mathematical}，
\begin{equation}
P_n(x) = \frac{ (2n-1)xP_{n-1}(x) - (n-1)P_{n-2}(x) }{n},
\end{equation}
在上面的式子中，x, $P_n(x)$都是绝对值不超过1的数，只有n可能是一个很大的数，我担心遇到“两个很大的数相减”的问题，所以将上面的公式重写为
\begin{equation}
P_n(x) = xP_{n-1}(x) + (1-\frac{1}{n}) ( xP_{n-1}(x) - P_{n-2}(x) ),
\end{equation}
根据这个公式，可以编写程序，得到勒让德多项式的值。
对于我们常用的100阶以内的勒让德多项式，我相信这个迭代公式都是很安全的。

\subsection{勒让德多项式的根}
这里我们可以证明，勒让德多项式的根满足一个不等式，通过这个不等式，我们可以界定这些根的取值范围，从而用二分法找到它们\cite{zeros-of-legendre}。

首先证明一个定理，在此基础上，我们推出我们要的不等式。
这个定理是Sturm定理的一个推论。
\begin{shaded}
$f(x)$与$F(x)$是$a\leq x\leq b$上的连续函数，且$f(x)\leq F(x)$但$f(x)$不恒等于$F(x)$，$y(x),Y(x)$分别满足微分方程
\begin{equation}
y'' + f(x) y = 0, ~~ Y'' + F(x) Y = 0.
\end{equation}
另有
\begin{equation}
y(x)>0 ~~ in~~ a<x<b, y(a) = y(b) = 0,
\end{equation}
那么，$Y(x)$不可能满足：$Y(x)\geq 0$且不恒为零。
\end{shaded}
这个可以用反证法证明。
我们假设$Y(x)\geq 0$且不恒为零。

我们注意到
\begin{eqnarray}
( y'Y - y Y' )|^b_a = \int^b_a ( F(x) - f(x) ) Y(x) y(x) dx 
\end{eqnarray}
而上式右边必大于0，且$y(a) = y(b)=0$，所以有
\begin{eqnarray}
y'(b)Y(b) - y'(a)Y(a) > 0,
\end{eqnarray}
然而，因为$y(a)=y(b)=0$，$(a,b)$区间内$y(x)>0$，所以必有$y'(a)>0, y'(b)<0$。
$Y(b), Y(a) \geq 0$，所以上式是不可能成立的。

推出矛盾，所以命题得到了证明。

容易看出来，$Y(x)$也不可能满足：$Y(x) \leq 0$且不恒为零。

类似地，如果$y(a)=y(b)=0$，在$(a,b)$区间内$y(x)<0$，也可以证明，$Y(x)$不可能满足以下任意一个条件
\begin{itemize}
\item $Y(x) \geq 0$且不恒为零。
\item $Y(x) \leq 0$且不恒为零。
\end{itemize}

所以，{\bf Y(x) 在$(a,b)$内一定有至少一个零点。}

勒让德多项式$P_n(x)$满足微分方程
\begin{eqnarray}
(1-x^2)P''_n(x) - 2xP'_n(x)+n(n+1)P_n(x)=0,
\end{eqnarray}
如果写作$P_n(\cos \theta),0<\theta<\pi$，则$z=(\sin \theta)^{1/2}P_n(\cos \theta)$满足微分方程
\begin{eqnarray}
z'' + ((n+1/2)^2 + (2\sin \theta)^{-2})z = 0,
\end{eqnarray}
另外我们知道，$y = \sin (n+1/2)x$满足
\begin{equation}
y'' + (n+1/2)^2 y = 0,
\end{equation}
$y(x)$在$[0,\pi]$有零点
\begin{eqnarray}
0, \frac{\pi}{n+1/2}, \frac{2\pi}{n+1/2}, \cdots, \frac{n\pi}{n+1/2}.
\end{eqnarray}

利用上面的定理，我们知道，$z=(\sin \theta)^{1/2}P_n(\cos \theta)$在$n$个开区间
\begin{eqnarray}
0 < \theta < \frac{\pi}{n+1/2} < \theta < \frac{2\pi}{n+1/2} \cdots < \theta < \frac{n\pi}{n+1/2},
\end{eqnarray}
中都至少有一个零点。
而$(\sin \theta)^{1/2}$在这些开区间中都不为零，$P_n(\cos \theta)$在这些开区间中总共最多有$n$个零点（因为它是n阶多项式），所以，$P_n(\cos \theta)$确然有$n$个零点，在这些开区间中分别有一个，
\begin{eqnarray}
0 < \theta_1 < \frac{\pi}{n+1/2} < \theta_2 < \frac{2\pi}{n+1/2} < \cdots < \theta_n < \frac{n\pi}{n+1/2},
\end{eqnarray}

\subsection{求取高斯勒让德积分的零点与权重}
根据Numerical Recipes，高斯积分
\begin{eqnarray}
\int^1_{-1}f(x)dx \approx \sum^n_{i=1} \omega_i f(x_i),
\end{eqnarray}
其中$x_i$为$P_n(x)$在[-1,1]上的零点，权重$\omega_i$的公式为\footnote{Numerical Recipes in C （中文版）上(4.5.9)应该是错的，(4.5.16)即下面的公式是对的}\footnote{http://mathworld.wolfram.com/Legendre-GaussQuadrature.html 链接中有一点信息。}
\begin{eqnarray}
\omega_i = \frac{2}{(1-x_i^2)P'_n(x_i)^2}.
\end{eqnarray}
下面的程序基本来自Numerical Recipes，我做了一点修改。
程序从初始点$cos( \frac{i+0.75 }{n+0.5} \pi )$出发，用牛顿法寻找$P_n(x)$的根，其中用到了$P_n(x)$的几个递推式。
找到根以后，根据上面的公式求取权重。
\begin{lstlisting}
/*
 * Given [a,b], returns abscissas x[] and weights w[], based on zeros of P_n(x)
 * Modified from Numerical Recipes
 */
void GauLeg(double a, double b, int n, double *x, double *w){

        double eps=1E-14;
        int m,j,i;
        double z1,z,xm,xl,pp,p3,p2,p1;

        m=(n+1)/2;
        xm=0.5*(a+b);
        xl=0.5*(b-a);
        for(i=0;i<m;i++){
                z=cos(M_PI*(i+0.75)/(n+0.5));//starting guess
                do{                     //Newton's method to get the ith zero
                        p1=1.0;         // P_0
                        p2=0.0;         
                        for(j=0;j<n;j++){
                                p3=p2;  // P_{j-1}
                                p2=p1;  // P_{j}
                                p1= z*p2 + (z*p2-p3)*j/(j+1); 
                                // P_{j+1},based on recurrence relation:
                                // (j+1)P_{j+1} = (2j+1)P_j - j P_{j-1} 
                        }
                        pp=n*(z*p1-p2)/(z*z-1); // P'n(x_i), based on
                                    // (1-x^2)P'_n(x) = n P_{n-1}(x) - nxP_n (x)
                        z1=z;
                        z=z1-p1/pp; // Newton's method: s = x - f(x)/f'(x)
                }while(fabs(z-z1)>eps);
                x[i]=xm-xl*z;
                x[n-1-i]=xm+xl*z;
                w[i] = 2.0*xl/((1.0-z*z)*pp*pp); // formula of weights:
                                    // w[i] = 2/(1-x_j^2)/P'n(x_j)^2 for [-1,1]
                                    // multiplied by xl, for [a,b]
                w[n-1-i]=w[i];
        }
}
\end{lstlisting}

\section{应用实例：氘核}
氘原子核由1个质子和1个中子组成，我们可以将它们之间的相互作用简化为有限深势井，即
\begin{equation}
V(r_{pn}) = \left\{
\begin{aligned}
&-V_0, ~~ &r < a, \\
&0, ~~ &else,
\end{aligned}
\right.
\end{equation}
其中$r_{pn}$是质子中子之间的距离，$V_0$为MeV量级，$a$为fm量级。
两体系统可以由质心坐标、相对坐标来描述。
质心坐标描述了原子核的平动，比如氘原子在水分子中，随着水分子在空气中做热运动。
相对坐标则描述了原子核内部的结构：e.g.基态时质子中子的相对距离。

相对坐标的薛定谔方程可以写作
\begin{equation}
(-\frac{\hbar^2}{2\mu} \nabla^2 + V(r) ) \Psi = E \Psi,
\end{equation}
其中$\mu　\equiv \frac{ m_p m_n }{m_p + m_n}$为折合质量。

分离变量以后，可以得到径向方程
\begin{equation}
(-\frac{\hbar^2}{2\mu} \frac{ d^2 }{d r^2} + V(r) + \frac{\hbar^2}{2\mu} \frac{l(l+1)}{r^2} ) u(r) = E u(r),
\end{equation}
而波函数为$\Psi = \frac{u(r)}{r} Y_{lm}(\theta, \phi)$.

在我们关于$V(r)$的假定下，体系的基态为$l=0$的束缚态。$l=0$时，径向方程为
\begin{equation}
(-\frac{\hbar^2}{2\mu} \frac{ d^2 }{d r^2} + V(r)) u(r) = E u(r),
\end{equation}
在$r <a$时，上述方程为
\begin{equation}
\frac{ d^2 }{d r^2}u(r)　+ K^2 u(r) = 0,
\end{equation}
其中$K = \frac{ \sqrt{ 2 \mu (E+V_0) } }{\hbar}$。

在$r <a$时，
\begin{equation}
\frac{ d^2 }{d r^2}u(r)　+ K^2 u(r) = 0,
\end{equation}
其中$K = \frac{ \sqrt{ 2 \mu (E+V_0) } }{\hbar}$。
$r > a$时，$u(r)$满足方程
\begin{equation}
\frac{ d^2 }{d r^2}u(r)　- \kappa^2 u(r) = 0,
\end{equation}
其中$\kappa = \frac{ \sqrt{ - 2 \mu E } }{ \hbar}$。

所以
\begin{equation}
u(r) = \left\{　
\begin{aligned}
& A\sin (Kr) + B \cos (Kr), & 0<r\leq a, \\
& Ce^{-\kappa r }, & r>a
\end{aligned}
\right.
\end{equation}

由于$\Psi = \frac{u(r)}{r} Y_{lm}(\theta, \phi)$在$r\rightarrow0$时为有限值，在$r=a$处连续，且一阶导数连续，得到
\begin{eqnarray}
B &=& 0, \\
A\sin(Ka) &=& C e^{-\kappa a}, \\ 
AK\cos(Ka) &=& -C \kappa e^{-\kappa a}.
\end{eqnarray}
得到
\begin{equation}
K \cot(Ka) = - \kappa,
\end{equation}
即本证能量为下面的函数的零点，
\begin{equation}
f(E) = \sqrt{2\mu(E+V_0)} cot( \sqrt{2\mu a^2 (E+V_0) } / \hbar ) + \sqrt{-2\mu E}. 
\end{equation}

取自然单位制$\hbar = c = 1$, $fm \approx 5 GeV^{-1}$，折合质量为$\mu = 0.5GeV$，
\begin{equation}
f(E)　(MeV) = \sqrt{1000(E+V_0)} cot( \sqrt{0.025 a^2 (E+V_0) } ) + \sqrt{-1000 E}, 
\label{eqn:deuteron-f(E)}
\end{equation}
其中$a$的单位为$fm$，$E, V_0$的单位为$MeV$。

如果取$a=2fm$，$V_0=35MeV$，$f(E)$的零点对应的是氘核的结合能，约2.2MeV，如图\ref{fig:deuteron}。
\begin{figure}
\centering
\includegraphics[width=0.8\linewidth]{deuteron}
\caption{Deuteron s-wave bound state: $f(E)$ in (\ref{eqn:deuteron-f(E)}) v.s. E in MeV, with $V_0=35MeV, a=2fm$。}
\label{fig:deuteron}
\end{figure}

